

clear mata
use  "${replication}\data\for_census_analysis.dta", clear

global breps 0

global int "elecXroad"
global  infs "elec road elecXroad"

estimates clear

global controls_base ""
global absorb "i.s_code"
global cluster "sd_code"

global depvars "D_icw_irrigation_2010_2012 irr_11 any_fruitvegspice secc_inc_cultiv_share secc_asset_index secc_pov_rate_rural"

global controls  "pc01_vd_p_sch_d pc01_vd_m_sch_d pc01_vd_s_sch_d pc01_vd_s_s_sch_d pc01_vd_college_d ln_pc01_pca_no_hh ln_pc01_pca_tot_p pc01_pca_p_sc_share pc01_pca_p_st_share pc01_pca_p_lit_share pc01_vd_power_agr pc01_vd_power_dom pc01_vd_power_all pc01_vd_dirt_road pc01_vd_tar_road ln_pc01_vd_area  ln_tdist_10 ln_tdist_50 ln_tdist_100 ln_tdist_500 tot_irr_share_2001  h_cntr_1_2001_d ph_cntr_1_2001_d drnk_wat_f_1_2001_d post_off_1_2001_d phone_1_2001_d " 

do "${replication}\code\Z_cons_boot.do"

global tmp "${replication}\Temp\"

cap label var elec_by_2010 "Electricity"
cap label var road_by_2010 "Road"
cap label var elecXroad_by_2010 "Elec and road"


global inf1 elec
global inf2 road
global int elecXroad 
**** 
cap rm "$tmp/cons_${int}_$inf1.csv"
cap rm "$tmp/cons_${int}_$inf2.csv"
cap rm "$tmp/cons_${int}_${int}.csv"

gen delta_irr=irr_11-tot_irr_share_2001
		

local column=0
foreach var of varlist $depvars{
	global terms ""
	local column=`column'+1
	
	reghdfe `var'   elec_by_2010 road_by_2010 elecXroad_by_2010  , absorb( $absorb) cluster($cluster)
	
	qui estimate store noc_S`column'
	
	qui estadd ysumm, mean
	
	if ("`var'"=="ln_secc_cons_pc_rural" | "`var'"=="secc_pov_rate_rural" ) & ${breps}>0 {
		preserve
		
		keep shrid  `var' elec_by_2010 road_by_2010 elecXroad_by_2010  s_code sd_code
		if "`var'"=="secc_pov_rate_rural" {
			local base_var secc_pov_rate_
			
		}
		if "`var'"=="ln_secc_cons_pc_rural" {
			local base_var secc_cons_pc_
			
		}
		fmerge 1:1 shrid using  "${replication}\data\shrug_rural_cons_boot.dta" , keep(match) keepusing(`base_var'*) nogen
		if "`var'"=="ln_secc_cons_pc_rural" {
			foreach tvar of varlist `base_var'* {
				replace `tvar'=ln(`tvar')
			}
		}
	
		cap rm $tmp/cons_boot_`var'_${inf1}.csv
		cap rm $tmp/cons_boot_`var'_${inf2}.csv
		cap rm $tmp/cons_boot_`var'_${int}.csv
		
		
		cons_boot, outfile("$tmp/cons_boot_`var'") name("`base_var'") spec("reghdfe `base_var'BOOTSTRAPNUM    elec_by_2010 road_by_2010 elecXroad_by_2010  [fw = weight_BOOTSTRAPNUM], absorb( $absorb) cluster($cluster)")  num_bs(${breps})
		
		
		restore
		
		estimates table noc_S`column'
		
		estimates restore noc_S`column'
		matrix V_temp=e(V)
		matrix V_temp[1,1]=(${se_${inf1}})^2
		matrix V_temp[2,2]=(${se_${inf2}})^2
		matrix V_temp[3,3]=(${se_${int}})^2
		erepost	V=V_temp
		estimate store  noc_S`column'
		
	}
	
	estimate restore  noc_S`column'
	
	lincomest elec_by_2010+road_by_2010+elecXroad_by_2010
	estimates store comb_noc_`column'
	
}

estimates table _all , keep(*by*) se stats(N $terms) 

global controls_missing ""

cap drop *missing
if "$controls"!="" {
	foreach var of varlist $controls {
		disp "`var'"
		qui gen `var'_missing=(`var'==.)
		qui replace `var'=0 if `var'==.
		cap gen `var'_SQ=`var'*`var'
		global controls_missing "$controls_missing `var'_missing "
	}
}

local column=0
foreach var of varlist $depvars{
	global terms ""
	local column=`column'+1
	disp "doing pds on `var'"

	pdslasso `var'      $inf1_by_2010 $inf2_by_2010 ${int}_by_2010 ($controls $controls_missing $absorb) , partial($controls_missing $absorb) cluster(sd_code)
	global selected_controls=e(xselected)
	global selected_controls_`column'="$selected_controls"
	global selected_controls_missing ""
	foreach control in $selected_controls {
		global selected_controls_missing "$selected_controls_missing `control'_missing"
	}
	
	reghdfe `var'     elec_by_2010 road_by_2010 elecXroad_by_2010  $selected_controls $selected_controls_missing  , absorb( $absorb) cluster($cluster)
	
	qui estimate store pds_`column'
		
	qui estadd ysumm, mean
	
	if ("`var'"=="ln_secc_cons_pc_rural" | "`var'"=="secc_pov_rate_rural" ) & ${breps}>0 {
		
		preserve
		
		keep shrid prop_tow_code IP_code `var' $inf1_by_2010 $inf2_by_2010 ${int}_by_2010  $controls $selected_controls  $selected_controls_missing s_code sd_code
		if "`var'"=="secc_pov_rate_rural" {
			local base_var secc_pov_rate_
			
		}
		if "`var'"=="ln_secc_cons_pc_rural" {
			local base_var secc_cons_pc_
			
		}
		fmerge 1:1 shrid using  "${replication}\data\shrug_rural_cons_boot.dta" , keep(match) keepusing(`base_var'*) nogen
		if "`var'"=="ln_secc_cons_pc_rural" {
			foreach tvar of varlist `base_var'* {
				replace `tvar'=ln(`tvar')
			}
		}
		
					
		cap rm $tmp/consp_boot_`var'.csv

		
		cons_boot, outfile("$tmp/consp_boot_`var'") name("`base_var'") spec("reghdfe `base_var'BOOTSTRAPNUM    ${inf1}_by_2010 ${inf2}_by_2010 ${int}_by_2010  $selected_controls  $selected_controls_missing [fw = weight_BOOTSTRAPNUM], absorb( $absorb) cluster($cluster)")  num_bs(${breps})
		
		store_est_tpl_boot using "${tmp}/consp_${int}_${inf1}.csv", infile("$tmp/consp_boot_`var'_${inf1}.csv") name("`base_var'") 
		store_est_tpl_boot using "${tmp}/consp_${int}_${inf2}.csv", infile("$tmp/consp_boot_`var'_${inf2}.csv") name("`base_var'") 
		store_est_tpl_boot using "${tmp}/consp_${int}_${int}.csv", infile("$tmp/consp_boot_`var'_${int}.csv") name("`base_var'") 
		restore

		estimates restore pds_`column'
		matrix V_temp=e(V)
		matrix V_temp[1,1]=(${se_${inf1}})^2
		matrix V_temp[2,2]=(${se_${inf2}})^2
		matrix V_temp[3,3]=(${se_${int}})^2
		erepost	V=V_temp
		estimate store  pds_`column'
		
		
		

		
	}
	
	estimate restore  pds_`column'
	lincomest ${inf1}_by_2010+${inf2}_by_2010+${int}_by_2010
	estimates store comb_pds_`column'
}

estimates table noc_*  , keep(*by*) se stats(N $terms r2) 
estimates table pds_*  , keep(*by*) se  stats(N $terms r2)
estimates table pds_*  , keep(*by*) star(.1 .05 .01)  stats(N $terms r2)

estimates table comb_noc_*  ,  se stats(N $terms r2) 
estimates table comb_noc_*  , star(.1 .05 .01) stats(N $terms r2) 
estimates table comb_pds_*  ,  se  stats(N $terms r2)
estimates table comb_pds_*  ,  star(.1 .05 .01)  stats(N $terms r2)


